 @include "glee"

 function numberp(x) {
     return x ~ /^[+-]?([0-9]+[.]?[0-9]*|[.][0-9]+)([eE][+-]?[0-9]+)?$/
 }
 function resetSeed(    seed) {
     if (seed) return(srand(seed))
     if(Seed)  return(srand(Seed))
     return srand(systime() + PROCINFO["pid"]) 
 }
 function sum(a, n,i) {
     for(i in a) n += a[i]
     return n
 }
 function a2cdf(a,cdf,    some,n,i,j,sum,memo) {
     for(i in a) {
	 some  = rand() * Pinch
	 a[i] += Pinch
	 memo[a[i]] = i
     }
     cdf[0] = n = asort(a,sorted)
     sum = cdf[1] = sorted[1];
     for(i=2;i<=length(sorted);i++) { 
	 sum += sorted[i];
	 cdf[i] = cdf[i-1] + sorted[i];
     }
     for(i in cdf)  
	 cdf[i] = cdf[i] / sum
     for(i=1;i<=n;i++)
	 cdf[-1*i] = memo[sorted[i]]
 }
 function pick(cdf,   i) {
     i = pick1(cdf,rand(),cdf[0])
     return cdf[-1*i]
 }
 function pick1(cdf,u,n,  i) {
     for(i=n-1;i>=1;i--) if (u > cdf[i]) return i+1;
     return 1;
 }	 
 function norm(x,min,max) { 
     return (x - min)/ (max - min + Pinch) 
	 }
 function stddev(n,sum,sumSq,   tmp) {
     tmp = (sumSq-((sum*sum)/n))/(n-1)
     tmp = tmp < 0 ? 0 : tmp
     return sqrt(tmp)
 }
 function gaussianPdf(m,s,x) {
     s = s + Pinch
     return 1/sqrt(2*Pi*s^2)*Ee^(-1*(x-m)^2/(2*s^2))
 }
 function round(x)      { 
     return int(x<0 ? x-0.5 : x+0.5) 
 }
 function abs(n) { return n < 0 ? -1*n : n }
 function odd(n) { return n % 2 } 
 function even(n) { return (n % 2) + 1 } 
 function median(a,    sorted,i,a1,n,half) {
     if (!sorted)  { for(i in a) a[i] += 0; asort(a,a1); return median(a1,1) }
     n = length(a)
     half = int(n/2) 
	 return odd(n) ? a[half] : (a[half] + a[half+1])/2 
	 }
 function array(a) {
     split("",a)
 }
 function push(a,v,old) {
     old      = length(a)
     a[old+1] = v
     return v
 }
 function pop(a,  old,v) {
     old = length(a)
     v   = a[old]
     delete a[old]
     return v
 }
 function pretty(x) {
     gsub(_,",",x)
     return x
 }
 function trim(x) {
     sub(/^[ \t]*/,"",x)
     sub(/[ \t]*$/,"",x)
     return x
 }
 function s2a(str,a,  sep,blankp,  tmp,n,i) {
     sep = sep ? sep : ","
     blankp = blankp ? blankp : 0
     n = split(str,tmp,sep)
     for(i = 1;i<=n;i+=2) 
	 if (blankp)
	     a[tmp[i]] = tmp[i+1]
         else a[trim(tmp[i])] = trim(tmp[i+1]);
     return int(n/2)
 }
 function o(a, str,control, i,com) {
     str = str ? str : "array"
     com = control ? control : " -n  -k 2" 
     com = "sort " com  " #" rand(); # ensure com is unique
     for(i in a) print oprim(str,i,a) | com;
     close(com) 
 }
 function oprim(str,i,a,  j) {
     j=i
     gsub(SUBSEP,",",j) 
     return str "[ " j " ]\t=\t [ " a[i] " ]"
 }
 function p2a(pairs,a,   sep,blankp,     n,tmp,i) {
     sep = sep ? sep : "/,/"
     n=split(pairs,tmp,sep)
     if (!blankp) 
	 for(i in tmp) 
	     tmp[i]=trim(tmp[i])
     for(i=1;i<=n;i=i+2) 
	 a[tmp[i]] = tmp[i+1] 
 }
 function rank(data0,ranks,					\
               data,starter,n,old,start,skipping,sum,i,j,r) {
     starter="someCraZYsymBOL";
     for(i in data0) data0[i] += 0
     n     = asort(data0,data)  
     old   = starter
     start = 1;
     delete ranks;
     for(i=1;i<=n;i++) {
	 skipping = (old == starter) || (data[i] == old);
	 if (skipping) {
	     sum += i 
	 } else {
	     r = sum/(i - start)
	     for(j=start;j<i;j++) 
		 ranks[data[j]] = r;
	     start = i;
	     sum   = i;
	 }
	 old=data[i]
     }
     if (skipping) {
	 ranks[data[n]] = sum/(i - start)
     } else {
	 if (! (data[n] in ranks))
	     ranks[data[n]] = n #int(r+1.4999)
     }
 }
 function criticalValue(conf) {
     conf = conf ? conf  : 95
     if (conf==99) return 2.326
     if (conf==95) return 1.960 
     if (conf==90) return 1.645
 }
 function minMax(col,x,missing,min,max) {
     if(x == missing) return 1
     if(col in min) {if(x < min[col]) {min[col]=x}} else {min[col]=x}
     if(col in max) {if(x > max[col]) {max[col]=x}} else {max[col]=x}
 }
 function readHeader(col,x,       name,indx,num,sym,dep,indep,	\
		     ignorep,nump,depp) {
     ignorep    = gsub(/\?/,"",x) 
     nump       = gsub(/\$/,"",x)
     depp       = gsub(/\!/,"",x) 
     indx[x]    = col
     name[col]  = x
     if (!ignorep)  {
	 nump ? num[col] : sym[col]
	 depp ? dep[col] : indep[col]
     }}
 function discretize(col,rows,data,bins,want1,b1,c1,want2,b2,c2,	\
		     all,n,row,sorted) {
     for(row=1;row<=rows;row++)  
	 all[row] = data[row,col]
     n    = asort(all,sorted)
     if(want1) efd(col,sorted,n,bins,b1,c1)
     if(want2) ewd(col,sorted,n,bins,b2,c2)
 }
 function efd(col,sorted,n,bins,b,c,   row,all,bin,i,jump,last) {
     bin  = last = 0
     i = jump = int(n/bins)
     while(i < n) {
	 if (sorted[i] == sorted[i+1]) {
	     i++
	 } else {
	     bin++
	     b[col,bin] = sorted[i]
	     c[col,bin] = i - last
	     last = i
	     i += jump }}
     if ((n-last) > (c[col,bin]/2)) 
	 bin++
     b[col,bin] = sorted[n]
     c[col,bin] += n - last
 }
 function ewd(col,sorted,n,bins,b,c,       min,max,jump,i,j) {
     min  = sorted[1]
     max  = sorted[n]
     jump = (max - min)/bins + 0.00000001
     for(i=1;i<=bins;i++)
	 b[col,i] = min + round(i* jump)
     j=1
     for(i=1;i<=n;i++) {
	 if(sorted[i]> b[col,j]) j++
	 c[col, j]++
     }
 }
 function asBin(col,x,bins,breaks,  bin) {
     for(bin=1;bin<=bins;bin++)
	 if ((col,bin) in breaks)
	     if (x <= breaks[col,bin])
		 return bin
     return 1
 }
 function nchars(n,char,   out) {
     while(n> 0) {n--; out = out char}
     return out
 }
 function fyi(x) {
     print "#FYI: " x >>"/dev/stderr"
 }


